In [ ]:
# Bibliotecas
import os                               # para manipulação de arquivos
import numpy as np                      # para operações numéricas
import librosa                          # para processamento de áudio
import matplotlib.pyplot as plt         # para criar gráficos e visualizações
import IPython.display as ipd           # para reproduzir áudio diretamente no Jupyter Notebook
from scipy.signal import stft, istft    # short term fourier transform
import pandas as pd                     # para criar tabelas
import seaborn as sns

1. Carregamento e Pré-processamento do Áudio¶

In [15]:
# Constantes otimizadas para análise de áudio
taxa_amostragem = 22050 # 22050 Hz é um valor de referência
nperseg = 1024          # valor de referência

# Carregamento de um audio de exemplo
audio, sr = librosa.load("audio/Confutatis.wav", sr=taxa_amostragem)  

ipd.Audio(data=audio, rate=taxa_amostragem)
Out[15]:
Your browser does not support the audio element.

2. Cálculo do Espectrograma (STFT)¶

In [16]:
# Calcular STFT 
frequencias, tempos, matriz_stft = stft(audio, fs=taxa_amostragem, nperseg=nperseg) 
espectrograma_db = 20 * np.log10(np.abs(np.abs(matriz_stft)) + 1e-10) # conversão para decibels
fase = np.angle(matriz_stft) # necessário para a reconstrução do áudio

3. Decomposição SVD e Análise Espectral¶

In [17]:
# Função para desenhar espectrograma
def plotar_espectrograma(S_db, sr, title="", ax=None):    
    img = librosa.display.specshow(S_db,
                                   sr=sr,
                                   hop_length=512,
                                   x_axis='time',
                                   y_axis='log',                                                                         
                                   cmap='inferno',
                                   vmin=-80, vmax=0,
                                   ax=ax)
    cbar = plt.colorbar(img, ax=ax, format='%+2.0f dB')
    cbar.set_label('Magnitude (dB)')
    if ax:
        ax.set_title("Espectograma"+" {}".format(title))
        ax.set_ylabel('Frequência [Hz]')
        ax.set_xlabel('Tempo [s]')
    else:
        plt.title("Espectograma"+" {}".format(title))
        plt.ylabel('Frequência [Hz]')
        plt.xlabel('Tempo [s]')

# Criar subplots lado a lado
plt.figure(figsize=(10, 4))

# Plotar espectrograma original
plotar_espectrograma(espectrograma_db, taxa_amostragem)

# Ajustar layout
plt.tight_layout()
plt.show()
No description has been provided for this image
In [18]:
def decomposicaoSVD(espectrograma, k=None):
    U, S, Vt = np.linalg.svd(espectrograma, full_matrices=False)
    U = U[:, :k]
    S = np.diag(S[:k])
    Vt = Vt[:k, :]
    return U, S, Vt
In [64]:
# Decomposição SVD do espectrograma
_, matriz_val_singulares, _ = decomposicaoSVD(espectrograma=espectrograma_db)

# Extrair os valores singulares
valores_singulares = np.diag(matriz_val_singulares)

# Calcular estatísticas
media = np.mean(valores_singulares)
desvio = np.std(valores_singulares)
valor_max = np.max(valores_singulares)
valor_min = np.min(valores_singulares)
mediana = np.quantile(valores_singulares,0.5)
quartil_1 = np.quantile(valores_singulares,0.25)
In [69]:
# Decomposição SVD do espectrograma
_, matriz_val_singulares, _ = decomposicaoSVD(espectrograma=espectrograma_db)

# Extrair os valores singulares
valores_singulares = np.diag(matriz_val_singulares)

# Calcular estatísticas
media = np.mean(valores_singulares)
desvio = np.std(valores_singulares)
valor_max = np.max(valores_singulares)
valor_min = np.min(valores_singulares)
mediana = np.quantile(valores_singulares,0.5)
quartil_1 = np.quantile(valores_singulares,0.25)
In [220]:
# Cálculo da proporção acumulada
proporcao_acumulada = np.cumsum(valores_singulares) / np.sum(valores_singulares)

# Encontrar o menor k tal que a proporção acumulada ≥ 95%
k_95 = np.argmax(proporcao_acumulada >= 0.95)

# Plot
plt.figure(figsize=(10, 4))
plt.plot(proporcao_acumulada, color='navy', linewidth=2, label='Energia acumulada')

# Linha vertical no k_95
plt.xticks(list(plt.xticks()[0]) + [k_95])

# Linha horizontal em 95%
plt.axhline(0.95, color='red', linestyle='--', linewidth=1.2, label='95% da energia')

# Rótulos e ajustes
plt.xlabel('Número de valores singulares $k$')
plt.ylabel('Proporção acumulada normalizada')
plt.title('Proporção acumulada dos valores singulares')
plt.grid(True, linestyle='--', alpha=0.6)
plt.legend()
plt.tight_layout()
plt.show()
No description has been provided for this image
In [235]:
# Paleta com alto contraste e baixa saturação
cor_media = '#5e81ac'     # Azul acinzentado
cor_mediana = '#a3be8c'   # Verde oliva suave
cor_quartil = '#d08770'   # Laranja queimado
cor_min_max = '#b48ead'   # Roxo malva apagado

# Plot
plt.figure(figsize=(10, 4))
sns.histplot(valores_singulares, kde=True, log_scale=True, color='#d8dee9', edgecolor='gray')

# Linhas verticais com legenda
plt.axvline(media, color=cor_media, linestyle='--', linewidth=2, alpha=0.9, label=f'Média (μ = {media:.2f})')
plt.axvline(mediana, color=cor_mediana, linestyle='-', linewidth=2, alpha=0.9, label=f'Mediana = {mediana:.2f}')
plt.axvline(quartil_1, color=cor_quartil, linestyle=':', linewidth=2, alpha=0.9, label=f'1º Quartil = {quartil_1:.2f}')
plt.axvline(valor_max, color=cor_min_max, linestyle='-.', linewidth=1.5, alpha=0.4, label=f'Máx = {valor_max:.2f}')
plt.axvline(valor_min, color=cor_min_max, linestyle='-.', linewidth=1.5, alpha=0.4, label=f'Mín = {valor_min:.2f}')

# Configurações visuais
plt.xlabel("Valores Singulares")
plt.ylabel("Frequência (escala log)")
plt.title("Distribuição dos Valores Singulares com Estatísticas")
plt.grid(True, linestyle='--', alpha=0.4)
plt.legend(loc='upper right')
plt.tight_layout()
plt.show()
No description has been provided for this image
In [166]:
# Criar o gráfico
plt.figure(figsize=(12, 4))
plt.plot(valores_singulares,".", color='#457b9d', linewidth=2, label="Espectro de Valores Singulares")
plt.grid(True)

# Linhas verticais de truncamento
locais_truncamento = [10,25,50,100,150,200,300,400]
for linha in locais_truncamento:
    plt.axvline(x=linha, color='red', linestyle='--', linewidth=0.75)

# Configurar os ticks do eixo x
xticks_custom = np.unique(np.concatenate((np.arange(0, len(valores_singulares), 100), locais_truncamento, [len(valores_singulares)])))
yticks_custom = np.unique(np.concatenate((np.arange(0, valor_max, 10000), [valor_max])))
plt.xticks(xticks_custom, rotation=45)
plt.yticks(yticks_custom)

# Legenda e títulos
plt.yscale('log')
plt.legend(['Valores Singulares', 'Locais de Truncamento'], loc='upper right')
plt.title(f"Valores Singulares (n={len(valores_singulares)})")
plt.xlabel('Quantidade de Valores Singulares')
plt.tight_layout()
plt.show()
No description has been provided for this image
In [94]:
locais_truncamento = [1, 10, 50, 100, 400]
n = len(locais_truncamento) + 1  # +1 para o original

ncols = 2
nrows = (n + 1) // 2  # número de linhas necessário

fig, axs = plt.subplots(nrows, ncols, figsize=(14, 3 * nrows))
axs = axs.flatten()  # transforma em array 1D para indexação sequencial

# 1. Espectrograma original
plotar_espectrograma(espectrograma_db, taxa_amostragem, title="Original", ax=axs[0])
axs[0].set_title("Original (completa)")

# 2. Espectrogramas truncados por SVD
for i, k in enumerate(locais_truncamento):
    idx = i + 1  # o índice 0 já foi usado para o original
    U, S, Vt = decomposicaoSVD(espectrograma=espectrograma_db, k=k)
    espectrograma_k = U @ S @ Vt
    plotar_espectrograma(espectrograma_k, taxa_amostragem, title=f"Truncado (k = {k})", ax=axs[idx])
    axs[idx].set_title(f"Truncado (k = {k})")

# Limpar eixos extras se houver
for j in range(len(locais_truncamento) + 1, len(axs)):
    axs[j].axis('off')

plt.tight_layout()
plt.show()
No description has been provided for this image

4. Reconstrução do Espectrograma e do Áudio¶

In [231]:
# Número de valores singulares a serem mantidos na decomposição SVD.
k = np.where(valores_singulares == np.quantile(valores_singulares, 0.25))[0][0]
U, S, Vt = decomposicaoSVD(espectrograma_db, k) 
espectrograma_resultante_db = U @ S @ Vt

# Convertendo de dB para magnitude linear 
magnitude_resultante = 10 ** (espectrograma_resultante_db / 20)

# Reconstruindo a matriz STFT complexa 
matriz_stft_resultante = magnitude_resultante * np.exp(1j * fase)

# Inverso da Transformada de Fourier de Curto Termo
_, audio_recuperado = istft(matriz_stft_resultante, fs=taxa_amostragem, nperseg=1024)

ipd.Audio(audio_recuperado, rate=sr)
Out[231]:
Your browser does not support the audio element.

5. Compressão vs Qualidade — Análise Quantitativa¶

In [250]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# --- Parâmetros ---
valores_k = [5, 10, 25, 50, 100, 150, 200, 300, 400, 450, 500, 513]
caminho_audio = "audio/Confutatis.wav"

# --- Matriz original ---
espectrograma_db 
linhas, colunas = espectrograma_db.shape
tamanho_original_elementos = linhas * colunas
tamanho_original_kb = espectrograma_db.nbytes / 1024

# --- Função de erro MSE ---
def erro_mse(original, reconstruido):
    return np.mean((original - reconstruido) ** 2)

# --- Pré-cálculo da SVD completa ---
U, S, Vt = decomposicaoSVD(espectrograma_db)

# --- Inicialização de armazenamento ---
resultados = []

# --- Loop sobre os valores de k ---
for k in valores_k:
    # Truncamento da SVD
    U_k = U[:, :k]
    S_k = S[:k, :k]
    Vt_k = Vt[:k, :]
    
    # Reconstrução do espectrograma_db
    espectro_reconstruido = U_k @ S_k @ Vt_k

    # Erro de reconstrução
    erro = erro_mse(espectrograma_db, espectro_reconstruido)

    # Tamanho em elementos (teórico)
    tamanho_elementos = k * (linhas + colunas + 1)
    taxa_elementos = tamanho_elementos / tamanho_original_elementos

    # Tamanho real em KB (efetivo)
    tamanho_kb = sum([matriz.nbytes for matriz in [U_k, S_k, Vt_k]]) / 1024
    taxa_kb = tamanho_kb / tamanho_original_kb

    # Armazenar resultado
    resultados.append({
        'k': k,
        'Erro MSE': erro,
        'Tamanho Comprimido (KB)': tamanho_kb,
        'Taxa de Compressão (KB)': taxa_kb,
        'Taxa de Compressão (Elementos)': taxa_elementos
    })

# --- Criar DataFrame com resultados ---
df_resultados = pd.DataFrame(resultados)

# --- Exibir a tabela ordenada ---
display(df_resultados.sort_values(by='k'))

# --- Gráficos lado a lado ---
fig, axs = plt.subplots(1, 2, figsize=(15, 5))

# Gráfico 1: Erro MSE
axs[0].plot(df_resultados['k'], df_resultados['Erro MSE'], marker='o', color='#f4a261')
axs[0].set_title('Erro de Reconstrução (MSE) vs k', fontsize=12)
axs[0].set_xlabel('Número de Componentes k')
axs[0].set_ylabel('Erro Quadrático Médio (MSE)')
axs[0].grid(True, linestyle='--', alpha=0.5)
axs[0].set_xticks(valores_k)
axs[0].tick_params(axis='x', rotation=45)

# Gráfico 2: Taxa de Compressão (KB)
axs[1].plot(df_resultados['k'], df_resultados['Taxa de Compressão (KB)'], marker='o', color='#457b9d')
axs[1].axhline(1.0, linestyle='--', color='red', alpha=0.6, label='Tamanho Original')
axs[1].set_title('Taxa de Compressão Real (KB) vs k', fontsize=12)
axs[1].set_xlabel('Número de Componentes k')
axs[1].set_ylabel('Taxa de Compressão (comprimido / original)')
axs[1].grid(True, linestyle='--', alpha=0.5)
axs[1].legend()
axs[1].set_xticks(valores_k)
axs[1].tick_params(axis='x', rotation=45)

plt.suptitle('Análise de Compressão via SVD', fontsize=14, y=1.05)
plt.tight_layout()
plt.show()
k Erro MSE Tamanho Comprimido (KB) Taxa de Compressão (KB) Taxa de Compressão (Elementos)
0 5 4.45e+01 146.07 0.01 0.01
1 10 3.96e+01 292.34 0.02 0.02
2 25 3.29e+01 732.32 0.05 0.05
3 50 2.80e+01 1469.53 0.11 0.10
4 100 2.15e+01 2958.59 0.21 0.21
5 150 1.65e+01 4467.19 0.32 0.31
6 200 1.24e+01 5995.31 0.43 0.42
7 300 6.47e+00 9110.16 0.65 0.63
8 400 2.71e+00 12303.12 0.88 0.84
9 450 1.35e+00 13928.91 1.00 0.94
10 500 2.34e-01 15574.22 1.12 1.05
11 513 1.24e-09 16005.20 1.15 1.07
No description has been provided for this image